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TEST  PARTICLE  CALCULATION  OF  ELECTRIC  CURRENTS 
IN  MAGNETIC  FIELD-REVERSED  REGIONS 

I.  INTRODUCTION 

Many  naturally  occurring  systems  can  be  thought  of  as  boundary  regions  between  oppo¬ 
sitely  directed  magnetic  fields,  where  the  dominant  magnetic  component  of  the  magnetic 
field  reverses  sign  and  where  “current  sheets”  are  formed.  Such  field  reversal  regions  occur 
in  laboratory  as  well  as  astrophysical  plasmas,  and  important  physical  phenomena  such  as 
reconnection  can  take  place. 

In  space  and  astrophysical  plasmas,  the  charged  particles  are  usually  collisionless.  A  simple 
example  of  a  field  reversal  which  models  the  Earth’s  magnetotail  is  the  quasi-neutral  sheet 
shown  in  Figure  1.  The  x-component  of  the  magnetic  field,  Bx{z),  reverses  sign  at  z  =  0. 
The  particle  motion  in  this  system  is  nonadiabatic  [1]  and  has  been  shown  to  undergo  chaotic 
scattering  [2,3].  Another  example  of  a  field  reversed  geometry  is  given  in  Figure  2.  Here  there 
is  an  X-type  neutral  line  at  x  =  0  and  z  =  0,  where  the  magnetic  field  vanishes.  The  particle 
motion  near  an  X-line  is  also  chaotic  [4,5]. 

Because  of  the  chaotic  particle  motion,  the  properties  of  plasmas  in  regions  containing 
current  sheets  are  nontrivial  to  study.  An  approach  that  has  been  fruitful  is  to  use  test  particle 
simulations  in  which  the  charged  particle  motion  is  integrated  in  prescribed  fields  [6-13]. 
Although  such  simulations  do  not  contain  complete  physical  interactions,  the  results  have 
often  been  successful  in  explaining  observed  satellite  data  [7,10] . 

A  key  property  of  magnetic  field  reversals  is  the  electric  current  profile  which  determines 
the  magnetic  field,  the  magnetic  energy  distribution,  and  particle  distribution  function,  the 
last  representing  the  internal  and  free  energy  distribution  of  the  plasma.  Issues  concerning 
current  distributions  have  been  investigated  using  test-particle  simulations  [13-19,  25].  For 
example,  Burkhart  et  al  [14]  modeled  the  structure  of  the  “dissipation”  region  in  collisionless 
reconnection.  In  this  work,  the  current  was  calculated  by  following  a  distribution  of  particles 
through  a  prescribed  field.  Using  the  calculated  current,  a  new  magnetic  field  is  calculated. 
The  process  was  iterated  until  Ampere’s  law  is  satisfied.  In  some  recent  papers,  self-consistent 
equilibrium  current  sheets  have  been  computed  in  the  magnetotail  geometry  [15,17,18].  In 
these  works,  the  cross-tail  current  density  is  again  constructed  by  binning  the  contributions 
from  each  test  particle  using  a  rectangular  (Cartesian)  grid  and  is  iterated  until  self-consistency 
is  achieved.  The  particle  orbits  are  computed  directly  from  the  equation  of  motion. 

An  important  feature  of  the  particle  motion  in  and  near  field  reversals  is  that  a  given  particle 
orbit  consists  of  alternating  segments  of  magnetized  and  unmagnedzed  motion  [1].  A  field 
reversal  thus  can  be  divided  into  spatial  regions  according  to  the  type  of  particle  motion. 
Consider  a  particle  of  energy  H  =  mv'^ll.  It  is  nearly  unmagnetized  while  in  the  region 
^  dj  =  where  Pz  is  the  particles  gyroradius  in  the  z-component  of  the  magnetic 

field  and  L  is  the  scale  length  of  the  reversal  region  in  the  x-component  of  the  magnetic  field, 
and  it  is  magnetized  for  |z|  ^  dj  [17]. 

In  this  report,  we  examine  in  detail  the  process  of  numerically  computing  electric  currents 
in  systems  in  which  the  particle  motion  alternately  undergoes  magnetized  and  unmagnetized 
motion.  In  addition,  we  describe  an  algorithm  for  approximately  separating  the  magnetic  field 
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into  the  part  due  to  the  plasmas  diamagnetism,  47rM,  and  the  part  due  to  the  actual  particle 
current,  H. 


n.  Single  Particles 


In  test  particle  simulations,  particle  motion  is  computed  in  prescribed  fields.  The  resulting 
equilibrium  profiles  (i.e.,  electric  currents,  particle  density,  and  pressure  tensor)  can  then  be 
used  to  calculate  the  fields,  which  in  turn  may  be  used  as  a  new  prescribed  field  in  the  test 
particle  simulation.  This  process  may  be  iterated  to  self-consistency.  Hence,  a  fundamental 
step  is  to  compute  the  contribution  of  each  particle  to  the  equilibrium  profiles  in  a  given  field. 
At  any  given  time  the  distribution  which  represents  a  single  particle  is  formally  given  by 


U{r,y,t) 
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Here  L  is  a  characteristic  scale  length  for  the  system  and  Hq  is  a  characteristic  frequency. 
Typically,  computation  of  physical  quantities  is  earned  out  using  Cartesian  coordinates  where 
we  have 


The  instantaneous  particle  and  current  densities  are  found  by  evaluating  the  zeroth  and 
first  velocity  moments  of  f^,  i.e.  n(r,  r)  =  J  d\  f{T,\,t)  and  j(r,r)  =  ^n(r,  r)U(r,  t)  = 
qJdv\f{T^  V,  t).  Defining  the  normalized  variables  f  =  r/L  and  v  =  v/ these  become 


n(f,t)  =  -^S[f-f(0] 


(3o) 


OLIU  Q 

Instead  of  evaluatingthe  pressure  tensorp(r,r)  =  m/t/v(v-U)(v-U)  /(r,v,0,  we  calculate 
the  second  velocity  moment  of /.p,i.e.  q(r,r)  =  m/<ivvv/(r,v,  t),  since  this  does  not  require 
a  priori  knowledge  of  the  average  flow  velocity  U.  Evaluating  this  integral,  we  find 

q(f,  t)  =  ^v(t)v(0  6[t  -  r(0].  (3c) 

It  is  trivial  to  show  that  q  may  be  written  in  terms  of  p,  j,  and  n  as, 

q  =  p  +  (m/^^n)jj. 
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(With  this  relation  in  mind,  we  shall  henceforth  simply  refer  to  q  as  the  pressure  unless 
otherwise  stated.) 

In  a  standard  particle  in  cell  code  it  is  the  instantaneous  density  and  current  which  is 
distributed  onto  a  grid  for  each  particle  at  every  time  step  [20].  The  resulting  total  density 
and  current  is  then  used  to  calculate  updated  fields  which  are  used  to  “push”  the  particles  for 
the  next  time  step.  In  self-consistent  test  particle  simulations,  on  the  other  hand,  a  particle 
is  “pushed”  through  the  system  in  a  prescribed  field  from  its  point  of  injection  to  its  point  of 
escape.  After  “pushing”  all  of  the  particles  in  a  source  distribution  through  the  system,  the 
average  currents  and  densities  each  particle  contributes  while  in  the  system  are  added  together 
to  determine  the  total  current  and  density  which  are  then  used  to  update  the  fields.  Thus,  the 
test  particle  code  is  very  well  suited  for  determining  time  stationary  solutions  with  very  good 
spatial  resolution  using  a  relatively  small  number  of  particles.  By  its  very  nature,  however, 
it  can  only  yield  information  regarding  time  evolution  in  the  sense  of  a  series  of  quasistatic 
equilibria. 

The  contribution  of  a  single  particle  to  the  equilibrium  quantities  is  obtained  by  averaging 
the  values  of  those  quantities  along  the  trajectory  for  the  time  T  the  particle  is  in  the  system  i.  e. 


«(r)  =  ^ 

f  n(f ,  t)dt 

Jo 

(4a) 

Kr)  =  j. 

f 

lo 

(4b) 

q(r)  =  ^ 

f  q(T,t)dt. 

/n 

(4c) 

(Note  that  the  value  of  T  will  in  general  be  different  for  each  particle.)  If  we  assume  that  the 
particle  phase  space  information  (i.e.  position  and  velocity)  are  known  at  the  evenly  spaced 
time  intervals,  =  nAt,  where  0  <  n  <  and  NAt  =  7,  we  may  approximate  the  integral  in 
(4a)  as  a  finite  sum.  Hence  the  average  single  particle  density  becomes 


n(f)  =  pD(f)  =  ^ 


1  ^ 

,  n=0 


(5a) 


Here  Sf  f.^  is  the  Kroneker  delta  function,  f „  =  f (nAt)  and  D(f )  is  the  normalized  single 
particle  density.  In  a  similar  manner,  the  average  current  density  and  pressure  which  may  be 
attributed  to  a  single  particle  may  be  written  as 


j(f)  =  ^j(f)  =  ^ 


i5]v(nAr)%. 

»=0 


(56) 


i53v(nA<)v(nA()«r,t. 

n=0 


(5c) 
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Where  J(f),  and  Q(f)  are  the  normalized  single  particle  current  density  and  pressure  tensor 
profiles  respectively.  Using  (5)  it  is  a  trivial  matter  to  distribute  the  single  particle  density  and 
current  onto  a  grid.  In  essence,  by  using  the  above  technique,  we  are  relating  average  residence 
times  of  particles  in  a  region  to  equilibrium  profiles. 

As  a  concrete  example,  we  consider  particle  motion  in  the  modified  Hams  model  of  the 
Earth’s  magnetotail  B  =  Bq  tanh(z/L)  x  +  z  (Fig.  1).  Here  Bq  is  the  asymptotic  value  of 
the  field  due  to  the  current  sheet,  L  is  the  characteristic  scale  length  of  the  field  reversal,  is 
a  constant  magnetic  field  normal  to  the  current  sheet,  and  x  is  in  the  sunward  direction.  The 
particle  motion  in  this  geometry  is  well  documented.  (A  review  can  be  found  in  reference 
[17].)  Since  the  fields  in  the  modified  Harris  model  are  translationally  invariant  in  the  x  and  y 
directions,  we  need  only  consider  variations  in  the  z-direction  when  calculating  the  current  and 
density  profiles.  Similar  techniques  may  be  used  for  X-line  geometries  (Figure  2),  however, 
in  that  case  we  must  use  a  two-dimensional  grid  since  there  is  variation  in  both  the  x  and  z 
directions. 

Assume  that  we  wish  to  distribute  the  current  and  density  profiles  onto  a  one  dimensional 
grid  in  the  z  direction  with  M-|- 1  equally  spaced  points  between  ±Zm<w  If  oii  the  time  step, 
the  particle  falls  between  the  m"'  and  (m  + 1  )"■  grid  points,  we  use  a  linear  interpolation  scheme 
to  distribute  the  Kroneker  delta  weight  function  onto  the  grid,  i.e.  we  add  to  the  density  and 
current  at  each  of  these  points  the  amounts 


An(i,)  = 

(6«) 

(6b) 

(6c) 

(6d) 

Aq(z„)  =  v(nA(MnA()J5^ 

(6e) 

Aq(z„+i )  -  v(nAr)v(nAr) ' 

(6f) 

Here  we  have  defined  Az  =  TzmaxiM.  This  procedure  is  foUowed  so  long  as  the  particle 
position  falls  between  ±Zmax-  Of  course,  if  we  are  to  use  the  top  (bottom)  grid  point,  we  must 
use  a  guard  cell,  i.e.,  if  the  particle  lies  between  Zmax  and  Zmax  +  Az,  {—Zmwc  and  —Zmwc  —  Az) 
we  must  add  to  the  density,  current  and  pressure  at  the  top  (bottom)  point  amounts  equal  to 
(6a),(6c)  and  (6e)  ((6b),(6d)  and  (6f))  where  now  z„  =  zm  iZm+\  =  zi)-  Finally,  in  order  to 
insure  unit  probability  of  finding  the  particle  in  the  system,  all  quantities  are  normalized  such 
that  the  height  integrated  density  (i.e.  the  sum  over  the  D(z)  contributions  at  each  grid  point) 
is  equal  to  one. 
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m.  Particle  distributions 


Having  found  the  equilibrium  profiles  due  to  a  single  particle  moving  through  the  field 
reversal  region,  we  now  wish  to  weight  these  single  particles  in  such  a  fashion  as  to  model  a 
source  distribution  of  particles.  In  this  context  there  is  some  ambiguity  as  to  how  to  perform 
the  reconstruction  based  on  what  is  assumed  to  be  known.  Here  we  shall  discuss  two  possible 
normalization  options.  In  both  cases,  we  assume  that  each  particle  represents  a  unit  volume  of 
phase  space  and  we  launch  shells  of  constant  energy  and  reconstruct  the  source  distribution  by 
properly  weighting  the  shells.  The  difference  between  the  two  noimalization  schemes  arise  in 

how  we  weight  the  particles  on  a  given  energy  shell. 

In  the  first  scenario,  (which  is  the  more  general  of  the  two)  we  assume  that  we  know 
the  density  of  the  incoming  distribution  and  that  the  density  of  the  outgoing  distribution  is 
determined  by  the  particle  dynamics.  In  this  case,  we  calculate  the  density  a  given  particle 
gives  to  the  top  grid  cell  when  it  is  initially  approaching  the  midplane  (assuming  the  particle 
is  launched  fi-om  above  the  field  reversal  region)  and  normalize  the  density  for  that  particle 
such  that  this  contribution  is  unity.  Note  that  this  normalization  implies  that  if  a  particle  leaves 
from  the  bottom  of  the  current  calculation  region,  the  density  of  the  top  cell  will  be  unity.  If 
the  particle  leaves  from  the  top,  however,  the  density  of  the  top  cell  will  be  greater  than  unity. 
Once  we  have  calculated  the  density  and  current  for  each  of  the  particles  on  a  given  energy 
shell,  all  of  the  single  particle  densities  and  currents  are  added  together  and  divided  by  the  total 
number  of  particles  that  enter  into  the  density  calculation  region  for  that  shell.  Note  that  in  this 
scenario,  we  may  either  symmetrize  the  current  sheet  by  assuming  equal  particle  sources  above 
and  below  the  midplane,  or  we  may  launch  separate  distributions  from  each  source  region. 

In  the  second  technique,  we  assume  a  priori  that  the  particle  sources  are  symmetric  and 
that  we  only  know  the  total  density  which  a  particle  contributes  to  the  top  grid  cell  is  equal  to 
unity.  In  this  case  we  symmetrize  the  single  particle  density  and  current  before  normaUzing 
them.  The  single  particle  equilibrium  profiles  are  then  normalized  such  that  the  density  in  the 
top  grid  cell  is  unity.  (Note  that  this  now  includes  particles  that  are  entering  and  leaving  the 
system  in  the  normaUzation.)  The  single  particles  profiles  for  a  given  energy  shell  are  added 
together  and  divided  by  the  total  number  of  particles  which  enter  the  current  calculation  region 

for  that  shell. 

In  general  the  first  technique  is  preferable  since  it  reduces  the  number  of  assumptions  made 
in  the  analysis.  Computationally  it  is  much  more  difficult,  however,  since  we  must  keep  track 
of  the  particle  while  it  is  initially  approaching  the  midplane  and  make  a  determination  as  to 
when  to  cut  off  the  calculation  of  the  normalization  factor  (i.e.  how  many  counts  the  particle 

contributes  to  the  top  grid  cell.) 

For  either  technique,  once  we  have  calculated  the  current  contribution  for  each  shell,  we 
normalize  the  profiles  by  their  relative  phase  space  volumes.  This  is  done  so  that  we  do  not 
need  to  space  the  shells  equally  in  velocity.  For  example,  suppose  we  have  P  energy  levels 
between  H  =  H„i„  and  H  =  H^.  The  shell  is  then  taken  to  represent  the  phase  space 

volume  (t.,-v,.i)/w  whercH=  ^m«2ai.dp=  1,2, The  phase  space  volume 
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represented  by  the  first  shell  is  simply  vi /v^.  Note  that  in  this  weighting  scheme  the  sum  of 
all  of  the  shells  has  unit  weight.  Finally  the  shells  are  weighed  by  the  desired  model  velocity 
distribution  function,  (e.g.  Maxwellian,  Kappa,  etc.)  and  summed  together  to  produce  the 
total  density  and  current  for  the  assumed  particles  sources. 


IV.  Application  TO  THE  Magnetotail 

In  the  case  of  the  magnetotail,  we  are  primarily  interested  in  the  y-component  of  the  current, 
since  this  is  the  component  which  produces  the  x-component  of  the  magnetic  field,  i.e.  the 
field  reversal.  Using  the  calculated  current,  we  compute  an  updated  magnetic  field  which  we 
then  use  as  the  input  field.  For  the  one-dimensional  case  discussed  here,  the  magnetic  field  is 
calculated  by  numerically  integrating  the  calculated  y-current  profile  in  the  z-direction  with 
the  numerical  constant  found  by  requiring  the  magnetic  field  to  vanish  at  the  midplane.  (For  a 
two-dimensional  case  such  as  the  X-line  geometry,  it  is  simpler  to  solve  for  the  vector  pqtential 
using  a  two-dimensional  Poisson  solver.)  In  updating  the  magnetic  field,  we  have  employed 
an  overall  weighting  factor  on  the  current  and  density  so  as  to  insure  that  the  amplitude  of 
the  asymptotic  magnetic  field  remains  unchanged.  The  process  is  iterated  until  the  input  and 
calculated  field  profiles  have  converged. 

Recently  it  has  been  found  that  if  all  three  components  of  the  current  are  kept,  an  antisym¬ 
metric  y-component  of  the  magnetic  field  grows  until  its  peak  value  is  on  the  order  of  the 
x-component  [21-23].  This  field  is  believed  to  be  the  test  particle  code  (time  independent) 
analog  of  a  wave  which  has  been  shown  to  develop  in  the  current  sheet  using  time  dependent 
hybrid  simulations  [23].  For  the  purpose  of  our  discussion  here,  it  suffices  to  note  that  if  the 
source  distribution  of  particles  is  highly  field  aligned,  this  wave  has  a  fast  growth  rate  and 
saturates,  with  By  on  the  order  of  Bx-  For  more  isotropic  source  distributions,  it  was  found  that 
the  growth  rate  is  much  slower  (on  the  order  of  half  an  hour  for  magnetotail  parameters)  and 
saturated  with  a  much  smaller  amplitude  [23].  Thus,  it  only  makes  sense  to  use  the  test  particle 
fonnulation  for  calculating  equilibria  for  the  cases  in  which  we  have  a  nearly  isotropic  source 
of  particles.  (It  is  also  important  to  note  that  in  order  to  successfully  run  the  hybrid  simulation, 
one  needs  to  input  the  self-consistent  magnetic  field  structure  calculated  using  the  test  particle 
code.  If  this  is  not  done,  the  hybrid  simulation  starts  too  far  from  equilibrium  which  results  in 
too  much  free  energy  for  wave  growth  and  a  resulting  disruption  of  the  equilibrium.)  Hence¬ 
forth,  we  will  suppress  the  y-component  of  the  magnetic  field  and  use  only  source  distributions 
that  produce  small  By. 

To  compare  out  results  with  observed  quantities,  we  must  convert  our  simulation  results 
into  physically  meaningful  units.  Using  equations  (5a-c),  we  see  that  the  total  particle  density, 
current  density  and  pressure  may  be  written  as 


1=1 

i=l 

where  C,  is  the  weighting  applied  to  the  i''*  particle  density,  Np  is  the  number  of  particles 
used  and  77  is  an  overall  weighting  factor  used  to  ensure  a  constant  asymptotic  magnetic  field 
strength.  We  now  need  to  relate  77,  L  and  our  dimensionless  density,  current,  and  pressure  to 
measured  values  of  the  asymptotic  magnetic  field  and  density. 

From  Ampere’s  law,  the  asymptotic  value  of  the  magnetic  field  may  be  related  to  the  height 

integrated  current  through  the  relation 


2Bo  =  «)^  r°° 

^  J  —  00 


In  our  computational  algorithm,  we  wish  to  we  maintain  a  constant  magnitude  for  the  asymp¬ 
totic  magnetic  field  {Le,,  B{—oo)/Bo  =  —  1  and  B(+oo)/J5o  =  +1)-  and  hence  we  require 


r°°  juiM/L) = 2 

J — oo 


Substituting  (9)  into  (8)  we  find  that  the  amplitude  of  the  magnetic  is  given  by 

Bo  =  4)rxlo22^nT  (10) 

Lj 

or  equivalently,  since  Oo  is  proportional  to  Bq  we  find  that 

77/L  =  5.19  X  10^^ -^m-^  (11) 

where  n  is  the  ratio  of  the  ion  to  the  proton  mass  and  Z  is  the  charge  state  of  the  ion.  Thus  we 
have  a  relationship  between  the  simulation  quantities  77/L  and  physically  determined  quantities. 
(Typically  we  assume  that  the  ions  are  protons  so  that  ^i  =  Z=\.) 

From  equation  (7a)  we  know  that  the  asymptotic  density  (or  alternatively  the  density  in  the 

top  cell)  is  given  by  ^ 

where  Zmax  is  the  rnaTcimnm  value  of  the  normalized  position.  If  we  require  that  htotizmax)  =  ”0 
where  no  is  the  physical  particle  density  given  in  particles  per  cubic  centimeter,  we  find  that 

JL  =  _  (12) 

L?  Djof(^max) 
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where  D,o,(  w)  is  the  computed  density  in  the  top  grid  cell  (output  from  the  code).  Hence 
we  have  a  relationship  between  the  physically  measured  particle  density  and  the  ratio  of  the 
simulation  quantities  -nil}  and  D,ot{zmwc) .  From  equations  (11)  and  (12)  it  is  now  possible  to 
uniquely  determine  the  values  of  r?  and  L  separately  for  any  given  physical  situation. 

In  space  physics  itis  often  useful  to  give  distances  in  terms  of  earths  radii  {Re  =  631x10  m). 
To  do  this  we  write  L  =  (tRe,  and  hence  from  (1 1)  we  have 

nl(T  =  3.307  X  10^^:^ 

Zj 


and  from  (12)  we  have  that 

n  _  ^0  >>3 

0-3  D,ot{Znujx)  ^ 

These  equation  may  be  solved  for  cr  to  obtain 


0-  =  0.036 


fj.no 


Z^^tot{Zmax) 


(13) 


In  essence,  equation  (13)  gives  us  the  scale  factor  for  changing  from  our  simulation  dimensions 
to  the  physical  dimensions  for  a  given  particle  density  and  ion  species. 

Using  the  above  results  it  is  now  a  trivial  process  to  determine  the  numerical  coefficients  to 
convert  our  simulation  density,  current  and  pressure  profiles  into  physically  meaningful  data, 

=  0.125{Bo/(t)  nA/m^.  (14) 

and  2 

=  7.97  X  10~^^(Bo^)  dyne/cm^.  (15) 

L 

In  equations  (13)-(15)  Bq  is  measured  in  nT  and  no  is  measured  in  cm“^. 

An  additional  quantity  of  interest  is  the  normalized  temperature,  f,  of  the  distribution  which 
is  input  to  the  code.  This  is  related  to  the  physical  temperature,  T,  through  the  relation 


f=  ITlmno^L^. 


(16) 


Solving  for  the  physical  temperature  in  keV,  we  find 

T  =  3.88^0  Vf(z2/Ai)  keV.  ( 17) 


In  Figure  3  we  give  an  example  the  self-consistent  density,  current,  magnetic  field  and 
pressure  profiles  calculated  using  the  above  algorithm  where  we  have  assumed  that  the  particle 
source  is  a  kappa  distribution  (k  =  4.5)  with  a  bulk  drift  in  the  tailward  direction  of  10% 
of  the  thermal  velocity  and  a  density  of  0.5cm~^,  the  asymptotic  magnetic  field  is  20nT,  and 

BjBo  =  0.1. 
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To  verify  the  existence  of  a  self-consistent  equilibrium,  the  density,  current  and  pressure 
profiles  must  satisfy  global  pressure  balance.  In  general  this  requirement  is  given  by  the  first 
moment  of  the  Vlasov  equation  which  yields 

V-(P  +  mNUU+|;:I-™)=0  (18) 

where  we  are  using  the  actual  pressure  P  =  m  /  tiv  (v  -  U)(v  -  U)  /{r,  u),  and  as  usual 
ATJ  =  /  dv  V  /(r,  v)  and  V  =  J  dv  /(r,  v)  are  the  particle  flux  and  particle  density  and  I  is  the 

unit  matrix.  In  ID  with  only  djdz  0,  Eq.  (18)  becomes 


P„ - -  =  const 

Alt 


(19a) 


yz 


ByBz 

Air 


=  const 


(19/.) 


Pzz  + 


Stt 


=  const. 


(19c) 


We  have  also  used  the  continuity  equation  and  the  symmetry  of  the  system  to  deduce  that 
NVz  =  0.  In  Figure  4  we  plot  equations  (19  a-c)  for  the  case  depicted  in  Figure  3.  Note  that 

all  three  equations  are  well  satisfied. 


V.  Effectsof  Diamagnetism 

An  interesting  aspect  of  the  field  reversed  geometry  is  the  effects  of  diamagnetism.  Even 
though  the  density  in  the  current  sheet  typically  only  varies  by  a  few  percent  in  the  cases  that 
produce  small  By,  we  find  that  it  is  in  these  cases  that  the  plasma  diamagnetism  becomes  a 
significant  effect.  This  is  because  the  particle  executes  magnetized  motion  for  lz|  ^  (pzL)  '  , 
and  unmagnetized  motion  for  jzj  ^  When  the  particle  motion  changes  character,  it 

contributes  a  diamagnetic  current  in  that  region.  As  we  will  show,  this  diamagnetic  current 
depends  on  the  asymptotic  pitch  angle  of  the  particles  with  larger  pitch  angles  producing 
larger  diamagnetic  currents.  In  this  section  we  will  discuss  the  diamagnetic  current  due  to  an 
individual  particle  and  then  discuss  the  ramifications  on  distributions  of  particles. 

Regardless  of  the  intervening  motion  (i.e.  chaotic  or  not),  the  shift  in  the  y-component  of 

the  guiding  center  is  given  by 


where  =  Bz/Bq  is  the  ratio  of  the  z  and  x  magnetic  fields  in  the  asymptotic  region, 
fj  _  is  the  normalized  energy  of  the  particle,  flo  =  qB^jmc  is  the  cyclotron 

frequency  in  the  x-component  of  the  magnetic  field,  and  /3i  (fii)  is  the  initial  (final)  pitch  angle 
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of  the  particle  as  it  enters  (leaves)  the  system.  From  (20),  we  may  calculate  the  average  current 
density  supplied  by  a  particle  to  be 

(j)  = 

where  At  is  the  time  required  for  the  particle  to  enter  and  leave  the  current  calculation  region 
[24], 

In  Figure  5a  we  show  the  trajectory  of  a  typical  orbit  and  in  Figure  5b  we  show  the  associated 
current  profile  in  the  y-direction  between  £  =  ±5  as  calculated  using  (6a)-(6d).  A  positive 
charge  moves  in  the  direction  indicated  by  the  arrow.  The  particle  is  launched  from  £o  =  6  with 
the  energy  H  =  and  an  initial  pitch  angle  of  100°  (using  the  field  conventions 

of  Figure  1).  In  its  interaction  with  the  neutral  sheet,  here  taken  to  be  the  time  the  particle  is 
in  the  current  calculation  region  {\z\  <  5),  the  y-component  of  the  particle  guiding  center  in 
the  asymptotic  regions  is  shifted  by  an  amount  Ay^^  =  3.03L  in  a  time  482JIq  ^ .  This  implies 
an  average  y-current  density  between  I  =  ±5  of  (j)  =  6.28  x  10  ^  qLQo.  If  we  average 
the  numerically  calculated  current  density  profile  between  z  =  ±5,  however,  we  obtain  a 
net  negative  total  current  of  -5.5  x  IQ-^qlHo.  This  somewhat  surprising  result  of  a  net 
negative  current  density  is  caused  by  the  plasma’s  diamagnetism  and  as  we  shall  see  is  strongly 
dependent  on  the  as5unptotic  pitch  angle  of  the  particle. 

To  understand  the  difference  between  the  numerically  computed  average  current  density  and 
the  value  given  by  gAy^^,/ Ar  the  size  of  the  current  calculation  region  is  increased  to  £  <  ±8 
so  that  the  particle  is  in  the  calculation  region  through  out  its  entire  orbit.  We  have  calculated 
the  Jy  profile  using  equations  (6a)-(6d)  at  several  times  to  see  how  the  current  is  deposited 
onto  the  grid.  Figure  6  shows  the  result  at  three  times:  (a)  T  =  200flQ  ,  (b)  T  =  400^0  ,  and 
(c)  T  =  600flQ  ^ .  The  important  feature  to  note  in  these  figures  is  the  development  of  current 
peaks  in  regions  associated  not  with  any  shift  in  the  y-component  of  the  guiding  center,  but 
rather  with  the  injection  and  exit  points  of  the  particles.  If  these  peaks  are  included  in  the 
numerically  calculated  average  current,  the  net  magnitude  of  the  current  density  is  positive  and 
in  agreement  with  the  value  calculated  from  (j)  =  qAygd^t.  Since  the  total  average  current 
(including  the  peaks  at  the  entry  and  exit  points)  is  in  agreement  with  the  expected  value,  this 
implies  we  have  deposited  a  net  negative  current  in  the  midplane  which  exactly  cancels  the 
positive  current  at  the  entry  and  exit  points.  This  negative  current  in  the  field  reversal  region 
is  in  addition  to  the  positive  current  due  to  the  meandering  motion  of  the  ion. 

The  formation  of  the  positive  current  peaks  at  the  entry  and  exit  points  and  the  deposition 
of  negative  currents  near  the  midplane  may  be  intuitively  understood  if  we  examine  a  plot  of 
V  vs.  £  (Figure  7a).  If  we  suppose  that  we  are  putting  the  current  into  bins  of  width  L,  then  it 
is  clear  that  the  finite  gyroradius  effects  in  the  regions  where  the  guiding  center  is  well  defined 
(Le.  from  points  a  to  b  and  from  points  c  to  d)  produce  a  net  positive  current  in  the  top  and 
bottom  grid  cells  but  tend  to  leave  net  negative  current  near  the  midplane. 

Another  way  of  understanding  the  development  of  the  current  peaks  at  the  entry  and  exit 
sites  is  to  examine  the  case  where  goes  to  zero.  Here  the  magnetic  field  is  purely  in  the 
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x-direction  but  we  are  still  griding  the  current  as  a  function  of  z  (Figure  8).  Exmnining  Figure  8, 
we  see  that  for  grid  cells  above  the  particle  guiding  center  there  are  only  positive  contributions 
to  the  y-current  density,  whereas  below  the  guiding  center  there  are  only  negative  contributions 
to  the  y-current  density.  Thus,  we  would  expect  to  develop  a  positive  current  peak  above  the 
guiding  center  position  and  a  negative  current  peak  of  equal  magnitude  below  the  guiding 
center  position.  If  the  two  peaks  are  added  together,  we  obtain  the  physically  reasonable  result 
of  no  net  current  in  the  y-direction.  The  same  difficulties  do  not  apply  to  the  x-current  density 
since  it  has  equal  positive  and  negative  contributions  both  above  and  below  the  particle  guiding 
center.  (This  formation  of  equal  and  opposite  current  peaks  above  and  below  the  guiding  center 
is  nothing  more  than  the  one-dimensional  manifestation  of  the  particles  diamagnetism.)  If  we 
now  allow  for  finite  Bz,  the  same  general  effects  will  occur  in  the  asymptotic  regions  (i.e.  far 
from  the  midplane),  but  now  the  peaks  will  tend  to  separate  as  we  let  T  increase  since  the 
particles  are  now  moving  down  the  field  line.  This  effect  is  illustrated  in  Figure  6a  where  we 
have  two  equal  and  opposite  current  peaks  separated  by  a  region  of  zero  current.  At  the  time 
depicted  in  the  Figure  6a,  the  particle  was  always  in  the  asymptotic  region  where  the  magnetic 
field  is  essentially  constant. 

In  appendix  A,  we  show  that  the  net  positive  y-velocity  deposited  around  the  entry  and  exit 
point  (or  conversely  the  net  negative  y-velocity  deposited  near  the  midplane)  for  a  particle  may 
be  approximated  by 


where 


The  quantity  {vy)p  is  related  to  the  current  density  in  the  peak  through  the  relation  ( Vy)p  = 
T{j  )p  where  Tis  the  total  time  between  when  the  particle  is  injected  into,  and  when  it  leaves 
the  current  calculation  region.  Figures  9a-9c  are  plots  of  the  approximate  analytic  expression 
(solid  line)  for  {vy)p  as  a  function  of  the  pitch  angle  for  ^  =  1.0  and  (a)  =  0.1,  (b)  =  0.2 

and  (c)  b„  =  0.4.  The  points  (•)  are  numerically  evaluated  by  “pushing”  the  particle  through 
the  assumed  field  and  keeping  track  of  the  y-velocity  deposited  around  the  entry  point  of  the 
particle.  We  note  that  for  nearly  field  aligned  particles  the  value  of  {vy)p  is  small  thus  indicating 
that  there  is  only  a  small  diamagnetic  current.  For  pitch  angles  near  90°,  however,  there  is 
substantial  diamagnetic  current  in  all  cases.  The  agreement  between  the  analytical  and  the 
numerical  results  is  in  good  agreement  throughout  the  parameter  space.  Note  the  breaks  in  both 
the  theoretical  and  numerical  curves.  Except  for  a  phase  shift  which  is  easily  understood  in 
terms  of  the  approximations  made  in  the  derivation  of  the  theoretical  expression,  the  locations 
and  number  of  breaks  are  in  good  agreement  with  the  numerical  results  (c/.  Appendix  A). 
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It  is  interesting  to  estimate  the  critical  pitch  angle,  /3c.  at  which  the  negative  diamagnetic 
current  deposited  in  the  vicinity  of  the  midplane  is  equal  to  the  positive  current  due  to  the 
meandering  motion  of  the  particle  in  the  field  reversal  region  ( i.e.  the  pitch-angle  at  which 
the  net  current  supplied  by  the  ion  is  in  the  wrong  direction  to  aid  in  the  production  of  the 
field  reversal).  As  a  lowest  order  approximation,  we  assume  that  as  the  particle  crosses  the 
midplane  all  of  its  energy  is  in  the  y-direction  and  that  it  remains  in  the  vicinity  of  the  midplane 
for  one  half  of  a  cyclotron  orbit  in  the  z-component  of  the  magnetic  field  [24].  (Note  that  this 
is  actually  an  over  estimate  of  the  meandering  motion  current  since  not  all  of  the  energy  will 
be  in  y-direction.)  This  implies  that  the  net  y-velocity  deposited  near  the  midplane  due  to  the 
meandering  motion  may  be  approximated  as 

{Vy)„  K  irVlk/bn  (23) 


Equation  (23)  is  shown  as  the  horizontal  dashed  lines  in  Figures  9a-9c.  Furthermore,  we 
assume  that  the  incoming  and  outgoing  pitch-angles  are  equal,  a  condition  which  has  been 
shown  to  be  nearly  the  case  for  Speiser-type  orbits  [18].  Equating  (23)  with  two  times  (22)  (to 
account  for  both  the  incoming  and  outgoing  diamagnetic  contributions),  it  is  easily  shown  that 
fic  is  a  solution  to  the  equation 


8  hz  sinj^c 
Try^l 


fc=0  \ 


1  - 


Irkb^cotfic 

y/l  bp' 


(24) 


In  Figure  10,  we  show  a  plot  of  /3c  as  a  function  of  b^,  (•)  as  given  by  (24).  We  also  show  on 
the  same  figure  the  value  of  Pc  found  by  using  the  actual  value  of  {vy)p  (i.e.,  the  solid  dots  in 
Figures  9a-9c)  instead  of  the  approximate  value  (i.e.,  the  solid  lines  in  Figures  9a-9c).  There 
are  two  important  facets  to  this  figure;  the  first  is  that  for  small  values  of  b^,  the  theoretical 
and  numerical  values  of  Pc  are  in  good  agreement,  the  second  is  that  for  larger  b^.  Pc  becomes 
smaller.  This  indicates  that  for  an  ensemble  of  particles  with  random  pitch  angles  the  ratio  of 
the  meandering  motion  current  to  the  diamagnetic  will  be  closer  to  unity  for  larger  values  of 

bz. 

An  important  conclusion  which  may  be  reached  from  the  above  discussion  is  that  if  the 
average  pitch  angle  of  a  source  of  particles  is  greater  than  Pc  the  diamagnetic  current  will  be 
larger  than  the  current  due  to  the  meandering  motion.  In  such  cases,  since  the  net  current 
is  in  the  wrong  direction  to  create  the  assumed  magnetic  field  reversal,  no  self-consistent 
solutions  exist.  In  essence  the  assumption  of  the  self-consistent  solution  to  the  magnetic  field 
and  the  chosen  distribution  function  are  mutually  exclusive.  If  on  the  other  hand  the  average 
distribution  function  has  a  large  drift  velocity  along  the  field  {vd  >  va),  the  average  pitch 
angle  is  much  less  the  Pc  and  the  diamagnetic  contribution  to  the  current  is  negligible.  In  this 
case  it  is  a  simple  matter  to  iterate  to  self-consistent  solutions  provided  that  By  is  negligible. 
[15,18,19,25].  As  mentioned  above,  one  must  be  careful  in  this  regime  since  By  is  typically 
non-n4ligible  and  grows  on  a  short  time  scale.  In  the  regime  reminiscent  of  the  quiet-time 
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magnetotail  where  the  field  aligned  drift  is  small  and  the  source  distribution  is  nearly  isotropic, 
the  average  pitch-angle  may  be  close  to  the  critical  pitch-angle.  Here  the  diamagnetic  current 
and  the  meandering  motion  current  are  of  the  same  order  magnitude  and  small  changes  in 
the  distribution  function  can  result  in  large  changes  in  the  self-consistent  solutions.  Previous 
attempts  to  calculate  the  current  sheet  structure  in  the  low  drift  velocity  regime  [15]  found  that 
if  the  drift  velocity  is  made  lower  than  some  critical  velocity,  then  no  equilibrium  exist.  It  was 
also  found  that  larger  values  of  bz  resulted  in  a  larger  value  of  the  critical  drift  velocity.  This 
loss  of  equilibrium  was  referred  to  as  the  current  sheet  catastrophe.  We  believe  that  what  the 
authors  were  actually  observing  was  the  transition  from  the  regime  in  which  the  meandering 
current  is  larger  than  the  diamagnetic  current  to  the  regime  in  which  the  opposite  is  true. 

VI.  Numerical  Algorithm 

In  this  section  we  describe  the  numerical  algorithm  used  to  separate  the  magnetic  field  into 
the  part  due  to  the  plasma  diamagnetism  and  the  part  due  to  the  actual  particle  current.  The 
simplest  possible  technique  is  to  follow  the  particle  guiding  center  and  to  evaluate  the  guiding 
center  drift  and  the  instantaneous  value  of  the  magnetic  moment  of  the  particle,  n  =  mv\/2B. 
Here  v±  refers  to  the  component  of  the  velocity  which  is  perpendicular  to  the  magnetic  field, 
B  is  the  local  magnitude  of  the  magnetic  field  and  the  direction  of  fi  is  antiparallel  to  the 
local  magnetic  field.  These  quantities  are  then  used  to  calculate  the  particle  current  and  the 
magnetization  vector  M  of  the  plasma  respectively.  (Note  that  the  magnetization  vector  M  is 
calculated  in  same  way  as  our  other  equilibrium  quantities,  i.e.  we  distribute  the  instantaneous 
magnetic  moment  of  each  particle  on  a  grid  and  then  add  them  together  with  the  appropriate 
weightings.)  This  technique  works  well  in  the  asymptotic  regions  where  the  guiding  center 
is  well  defined  but  fails  near  the  midplane  where  the  guiding  center  is  only  poorly  defined. 
We  therefore  use  a  hybrid  model  for  calculating  the  equilibrium  profiles.  In  the  regions  where 
the  guiding  center  is  well  defined,  we  transform  to  the  guiding  center  coordinates  (position 
and  velocity)  for  distributing  the  current  density,  density  and  magnetic  moment  onto  the  grid 
whereas  in  regions  where  the  guiding  center  is  not  well  defined  we  use  the  actual  particle 
position  and  velocity  for  distributing  the  current  and  density  onto  the  grid.  (See  appendix  B  for 
the  calculation  of  the  guiding  center  location.)  This  method  requires  that  we  follow  the  actual 
particle  motion  in  both  regions,  so  that  proper  phase  information  is  retained  for  the  particle  in 
changing  from  the  magnetized  to  the  unmagnetized  regions  and  vice  versa.  For  the  case  of  the 
magnetotail,  the  separation  between  the  “magnetized”  and  “unmagnetized”  regions  occurs  at 
the  separatrix  where  the  particle  changes  from  non-midplane  crossing  to  midplane  crossing. 
In  appendix  C,  we  show  that  the  phase  space  location  of  the  separatrix  is  approximated  by 

u/  <  [Py/M- {q/Mc)Ay{x,z)]‘^  -  [Py/M- iqlMc)Ay{x,z  =  0)]^.  (25) 

where  Ay  is  the  vector  potential  for  the  magnetic  field  and  Py  is  the  y  component  of  canonical 
momentum  for  the  particle.  The  condition  given  by  (24)  is  simply  a  statement  that  for  a  particle 
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to  be  classified  as  a  midplane  crossing,  the  energy  in  the  z-component  of  the  velocity  must  be 
larger  than  the  effective  potential  for  the  z-motion  evaluated  at  the  midplane. 

In  Figure  11,  we  show  the  density,  current  and  magnetic  field  profiles  calculated  using 
the  algorithms  described  in  this  paper  for  the  case  of  a  Maxwellian  with  a  bulk  drift  in  the 
tailward  direction  of  10%  of  the  ion  thermal  velocity  and  a  high  energy  isotropic  tail;  which 
comprises  10%  of  the  total  density.  Note  that  in  this  case,  (which  has  a  nearly  isotropic  source 
distribution)  the  diamagnetization  accounts  is  approximately  80%  of  the  magnetic  field  due 
to  the  particle  currents.  In  Figure  12  we  show  the  same  plots  but  using  a  Maxwellian  source 
distribution  which  has  a  tailward  drift  of  twice  the  ion  thermal  velocity  (i.e.,  a  highly  field 
aligned  distribution).  Note  that  in  this  case  the  plasma  diamagnetism  is  only  a  small  correction 
to  the  total  magnetic  field.  One  must  be  careful  in  interpreting  the  results  in  this  case;  however, 
since  we  would  expect  there  to  be  a  rapidly  growing  By  field  associated  with  the  highly  field 
aligned  distribution. 


vn.  Conclusions 

In  this  report  we  have  described  a  technique  for  calculating  density,  electric  current  and 
pressure  in  a  magnetic  field  reversed  region  using  a  test  particle  simulation.  The  essence  of 
the  technique  is  to  push  a  distribution  of  particles  through  a  model  magnetic  field  which  is 
fixed  in  both  time  and  space.  For  each  particle  the  density  current  and  pressure  are  laid  down 
a  on  grid.  After  all  of  the  particle  trajectories  have  been  calculated  we  then  calculate  the 
total  density,  current,  pressure  and  magnetic  field  profiles  that  the  assumed  distribution  would 
produce.  We  then  compare  the  assumed  and  the  calculated  magnetic  fields.  If  the  difference 
between  the  fields  is  sufficiently  small  the  simulation  is  stopped  and  we  check  to  make  sure 
that  global  pressure  balance  is  obtained  [Eq.  (19)].  Otherwise,  we  mix  the  assumed  and  the 
calculated  fields  and  use  the  combined  field  in  a  new  iteration  of  the  test  particle  code.  This 
technique  has  proven  to  be  very  effective  for  calculating  self-consistent  equilibrium  structures 
in  complex  systems  such  as  the  magnetotail  because  it  gives  very  good  spatial  resolution 
using  a  relatively  small  number  of  particles.  As  discussed  above,  however,  it  can  only  yield 
information  regarding  time  evolution  in  the  sense  of  a  series  of  quasistatic  equilibria. 

Using  this  method  we  have  calculated  self-consistent  solutions  for  the  equilibrium  structure 
of  the  Earth’s  magnetotail  during  quiet-times  which  are  in  good  agreement  with  structures 
observed  in  satellite  data.  An  interesting  aspect  of  the  solutions  is  the  effect  of  the  plasma 
diamagnetism.  Using  a  post-processing  algorithm  (Sections  V  and  VI),  we  have  separated  the 
magnetic  field  into  the  components  due  to  the  plasma  diamagnetism,  47rM,  and  the  part  due  to 
the  actual  particle  current,  H.  We  find  that  if  the  average  pitch-angle  of  the  source  distribution 
is  too  large,  the  magnitude  of  the  diamagnetic  current  is  large  than  the  magnitude  of  the  particle 
current  and  in  the  wrong  direction  to  produce  the  assumed  field  reversal.  Obviously,  in  such 
cases  the  two  assumptions,  ue.,  the  assumed  form  of  the  distribution  function  and  the  assumed 
magnetic  field  structure,  are  mutually  exclusive  and  no  self-consistent  solutions  exist.  We 
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further  showed  that  for  larger  values  of  the  average  pitch  angle  required  for  the  existence 
of  self-consistent  solutions  becomes  smaller,  i.e.  we  require  a  more  field  aligned  distribution 
function.  We  believe  that  these  results  explain  the  reported  current  sheet  catastrophe  in 
Reference  [15]. 
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APPENDIX  A 


An  interesting  feature  of  the  results  presented  in  this  paper  is  that  as  the  particle  distribution 
function  becomes  more  isotropic  Gcss  field  aligned),  the  diamagnetic  contribution  to  the  total 
cross-tail  current  becomes  larger  and  that  in  extreme  cases  it  is  large  enough  to  cancel  the 
cross-tail  current  due  to  the  ’’meandering  motion”  of  the  ions  as  they  cross  the  midplane. 
Obviously,  in  these  extreme  cases,  the  assumption  of  a  magnetic  field  reversed  topology  and 
the  assumed  distribution  function  are  mutually  exclusive,  since  the  calculated  current  is  either 
zero  or  in  the  wrong  direction  to  produce  the  desired  field  reversal.  It  is  also  observed  that 
a  larger  field  aligned  anisotropy  is  required  for  convergence  for  larger  values  of  than  for 
smaller  values  of  bz.  It  is  the  aim  of  this  appendix  to  given  a  theoretical  estimate  of  the  effects 
of  pitch  angle  and  b^  on  the  diamagnetic  current. 

Since  the  net  diamagnetic  current  of  a  particle  in  a  closed  system  must  be  zero,  we  may 
calculate  the  net  diamagnetic  current  that  any  given  particle  deposits  in  the  vicinity  of  the 
midplane  by  examining  the  diamagnetic  current  that  a  particle  deposits  as  it  enters  and  leaves 
the  system.  In  practice,  rather  than  calculating  the  actual  cross-tail  current  density,  we  calculate 
the  net  y-velocity  an  ion  deposits  as  it  passes  through  the  calculation  region.  Whereas  these 
only  differ  by  a  normalization  factor,  it  greatly  simplifies  the  analysis. 

To  determine  the  net  positive  y- velocity,  {vy),  a  particle  deposits  around  the  point  where 
it  is  launched  into  the  system  and  the  point  where  it  escapes  from  the  system  we  need  to 
calculate  the  z-position,  the  z-velocity  and  the  y- velocity  as  a  function  of  time.  In  terms  of  our 
normalized  units,  z  =  zjL,  =  Vz/{^L),  and  Vy  =  Vy/ {^IqL)  these  are  given  by 


z  =  zo  + 


bzTCOS  0  -H 


\/l  + 


;  sin  yS  sin  ( r ) 


(Al) 


2H 


bzCOsP  +  sin/3  cos( 1  +  bz^  r) 


Vy  =  \F2il  sin/?  sin 


(A2) 

(A3) 


where  H  =  (mv^/2) /(mflo  V)  is  the  normalized  energy,  r  =  r/i  is  the  normalized  position, 
^  is  the  pitch  angle  of  the  particle,  and  =  B„/ Bq  is  the  ratio  of  the  constant  z-component  of 
the  magnetic  field  with  the  asymptotic  value  of  the  x-component  of  the  magnetic  field. 

During  a  single  excursion  of  the  particle  above  the  z  =  zo  pl3ii®»  we  add  to  the  net  y- velocity 

in  this  region,  {vy)p,  an  amount 


(vyk)  =  f  sinP  sm{yf[  +  ^T)  dr  (A4) 

Jr2k 

where  T2k,  and  r2*+i ,  are  the  times  at  which  the  particle  cross  into  and  out  of  the  region  z  >  zo 
respectively.  In  actuality,  this  implies  that  T2k  and  t^jc+i  should  be  two  consecutive  solutions 
to  the  nonlinear  equation 
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b^TCOS p  +  .  ^  --  sin  /?  sin (^1"+^ r )  =  0  (A5) 

Vl+^z^ 

Rather  than  solve  (A5)  for  each  excursion  of  the  particle  above  the  z  =  zo  plane,  we  hold 
the  position  of  the  guiding  center  constant  in  performing  the  integration.  We  then  shift  the 
location  of  the  guiding  center  down  by  an  amount,  Az,  equal  to  what  it  would  drift  in  a 
complete  cyclotron  period.  This  allow  us  to  obtain  approximate  analytic  values  for  the  time 
limits  in  the  integration.  The  shift  in  the  z  location  between  two  successive  integration  is  given 

by 

^  x/ii cos/3  (A6) 

Note  that  if  >  Rg,  where  is  an  integer  and  Rg  is  the  normalized  gyroradius  of  the 
particle,  the  particle  will  not  go  above  the  z  =  zo  plane  and  will  thus  no  longer  contribute  to 
the  net  y-velocity  deposited  in  that  region.  Referring  to  Figure  (13),  we  see  that  in  a  given 
excursion  above  the  z  =  zo.  the  particles  phase  angle  goes  from  to  tt  -  Ok  and  hence  the 
time  integration  runs  from 

r2k  =  {l  +  (^7a) 

to 

T'ik+l  =  (1  +  (27r(^  +  1/2)  -  h) 

where 

dk  =  sin"^  r  2irkbz  J 

1-VT+^  J 

Carrying  out  the  integration  with  the  approximate  limits,  we  obtain 

(vyk)  =  Sin/3  cosSk  (A9) 


sin/3 

The  total  y-velocity  deposited  above  the  z  =  zo  plane  is  found  by  summing  the  y-velocity 
deposited  in  that  region  during  the  individual  excursions,  i.e., 

{Vy)tot  =  '^{^yk)- 
k=\ 
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Here  K  is  the  number  of  excursions  that  the  particle  makes  above  the  z  =  zo  plane  and  is  given 
by 


^:=Int 


^^^tan^ 


ItT&z 


(A12) 


Combining  equations  (AlO),  (All),  and  (A12)  yields  an  expression  for  the  net  positive  y- 
velocity  deposited  above  the  z  —  zo  plane.  In  addition,  we  have  an  equal  amount  of  positive 
y- velocity  deposited  symmetrically  below  the  z  =  ^  plane.  This  is  because  until  a  particle  has 
drifted  at  least  one  gyroradius  in  the  z-direction  the  negative  and  positive  contributions  to  the 
y-current  cannot  average  to  zero.  Thus,  we  see  that  the  total  y-velocity  in  the  peaks  located 
around  the  entry  and  exit  points  of  the  particle  is  given  by 


{Vy)p  =  2  {Vy)tOf  (^13) 

In  Figure  9  we  show  the  total  y-velocity  in  the  current  peaks  as  a  function  of  the  pitch  angle  as 
given  by  given  by  (A  13)  (solid  line)  and  as  determined  by  numerically  “pushing  the  particle” 
through  the  field  (•).  As  can  be  seen,  the  agreement  is  quite  good.  One  notes,  however, 
that  the  location  of  the  “breaks”  in  the  curves,  which  occur  whenever  the  particle  is  able  to 
add  an  additional  excursion  above  the  z  =  ^  plane,  fall  in  slightly  different  locations  for  the 
theoretical  and  numerically  determined  results.  For  the  theoretical  curves,  the  breaks  in  the 
curves  occur  at  those  pitch-angles  for  which  the  particle  is  able  to  make  an  additional  excursion 
above  the  z  =  zo  plane,  i.e.  when 

In  actuality,  the  breaks  in  the  curve  occur  at  those  pitch  angles  for  which  (A5)  picks  up  two 
additional  roots.  This  happens  whenever  the  amplitude  of  the  sinusoidal  variation  in  r  becomes 
sufficiently  large  that  an  additional  peak  intersects  the  linear  term  in  (A5).  Again,  the  actual 
values  require  a  solution  to  the  full  nonlinear  equation.  For  k=0,  this  is  easily  done  yielding 

ySo  =  tan“^(i)z).  (A15) 

It  is  important  to  note  that  for  predominantly  field  aligned  orbits,  i.e.,  those  with  pitch  angles 
less  than  A)  as  given  by  (A15)  there  are  no  solutions  to  (A5)  thus  indicating  that  there  is  no  net 
y-velocity  deposited  above  the  z  =  zo  pl£ine.  Note  that  for  pitch-angles  less  than  /3d  as  given 
by  (A  15)  the  numerically  determined  points  in  Figures  9a-9c  are  all  zero.  For  the  other  roots 
the  solution  is  not  as  simple.  As  a  lowest  order  approximation  to  the  actual  solutions  we  may 
approximate  the  roots  as  occurring  whenever  the  amplitude  of  the  sinusoidal  term  in  (A5)  is 
equal  to  the  value  of  the  linear  term  at  the  points  -v/I  +  ^r  =  (4/c  +l)Tr/2.  This  yields  the 
values  of  given  by 
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(A16) 


^*  =  tan  ^  Trfcj]. 

In  actuality,  the  solutions  occur  for  values  of  t  less  than  this  and  thus  the  values  of  fik  are 
slightly  off.  Referring  to  Figure  14,  however,  where  we  have  plotted  the  linear  term  in  (A5) 
along  with  three  examples  of  the  sinusoidal  term  with  amplitudes  given  by  (a)  Pq  as  given  by 
(A15)  and  (b)/?i  and  (c)  ^2  as  given  by  (A16)  we  see  that  (A16)  gives  an  excellent  approximation 
to  the  actual  value  of  Comparing  equations  (A14)  and  (A16)  we  see  that  there  is  a  small 
positive  shift  in  the  actual  location  of  the  Pk  over  the  approximate  location. 

APPENDIX  B 


In  this  appendix  we  explicitly  write  out  the  transformation  between  the  actual  particle 
position  and  the  guiding  center  position  for  magnetic  fields  of  the  form  B  =  Bq  /(z)x  + 

In  general,  the  transformation  is  given  by 

«  vxB 

{g/Mc)B^'  ^  ^ 

whereB  =  [Bof{z)]^+Bz^.  Using (Bl),itistrivialtoshowthatthe guiding centercoordinates 
Xgc,  Ygc,  and  Zgc  are  given  by 


Xgc  —  ^  4* 
Ygc  =  y  + 
Zgc  -z  + 


bzVy 

^o[f{z)'^  +  bz^] 

(B2a) 

f(z)  Vz  -1-  bz  Vx 

f{z)v, 

^o[f{z)^  + 

(Blc) 

Thus,  we  see  that  given  the  local  value  of  the  magnetic  field  and  the  instantaneous  velocity, 
it  is  a  simple  matter  to  calculate  the  instantaneous  guiding  center  position.  For  unmagne¬ 
tized  motion  in  the  field  reversal  region,  the  instantaneous  guiding  center  has  little  physical 
significance  or  practical  utility. 


APPENDIX  C 

In  sections  V  and  VI,  we  showed  that  in  order  to  determine  the  effects  of  plasma  diamag¬ 
netism  on  the  equilibrium  structure  of  the  magnetotail  we  needed  to  separate  the  orbit  into 
segments  which  are  midplane  crossing  and  segments  which  are  non-midplane  crossing.  In 
the  former,  when  calculating  equilibrium  profiles  we  distribute  the  actual  particle  positions 
and  velocities  onto  a  grid,  whereas  in  the  latter  we  distribute  the  guiding  center  positions  and 
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velocities  onto  a  grid.  In  this  appendix  we  present  an  approximate  expression  for  the  phase 
space  location  where  the  particle  changes  its  behavior  from  crossing  to  non-crossing. 

We  begin  by  considering  the  Hamiltonian  of  a  particle  in  the  modified  Harris  magnetic  field 

H{X,  Z,  V^,  V„  Py)  =  -H  U{X,  Z,  V;,,  Py)  (Cl ) 

where 

Uix,z,v,,Py)  =  ^Mv'^  +  ^MlPy/M- {q/Mc)Ay{x,z)]^.  (C2) 

is  the  effective  potential  for  the  z-motion,  Py  is  the  y-component  of  the  canonical  momentum 
and  Ay  is  the  vector  potential.  For  the  modified  Harris  magnetic  field,  Ay  may  be  written  as 

Ay(x,z)  =  — BoLln[cosh(z/L)] -I- BjX  (C3) 

If  hz  -C  1  the  oscillations  in  the  z-direction  are  much  faster  than  the  oscillations  in  the  x- 
direction.  As  a  first  approximation,  we  therefore  assume  that  during  a  single  z-oscillation  the 
values  of  both  and  x  are  constant.  This  approximation  is  not  strictly  valid.  Rather,  if  we  take 
Vx  to  be  a  constant  we  should  approximate  x  =  xq  +  f.  For  our  purposes,  however,  we  find 
that  the  simple  approximation  is  sufficient.  This  is  because  the  orbits  which  incur  significant 
error  using  the  standard  current  calculation  technique  have  large  pitch  angles  and  thus  small 
Vx  Setting  Vx  and  x  constant,  the  condition  that  a  particle  is  non  axis  crossing  is  given  by 

H{x,  z,  Vx,  Uz,  Py)  <  U{x,  z  =  0,  Vx,  Py)  (C4) 

i.e.,  the  energy  of  the  particle  must  be  grater  than  the  effective  potential  at  the  midplane.  This 
criterion  is  graphically  depicted  in  Figure  15  where  we  have  drawn  the  effective  potential  for 
{PyfM- n„x) /(fioB)  =  (a)  1,  (b)  0,  and  (c)  -1,  and  we  have  defined  =  qBjMc.  We  have 
neglected  the  Mv^jl  contribution  to  the  effective  potential,  since  it  only  provides  a  constant 
offset.  The  dashed  lines  correspond  to  orbits  of  different  energies  as  compared  to  effective 
potential  (c).  Particle  A  violates  the  inequality  in  (C4)  and  is  midplane  crossing.  Thus  we 
use  actual  particle  coordinates  for  calculating  the  equilibrium  profiles.  Particle  C  represents 
a  solution  to  the  inequality  (C4)  and  is  non-midplane  crossing.  Thus  we  would  use  guiding 
center  coordinates  for  calculating  the  equilibrium  profiles.  Particle  B  represents  the  equality 
condition  in  (C4)  and  is  at  the  transition  point  between  crossing  and  non-crossing  motion. 
Solving  (C4)  for  v^  we  find  that  the  condition  a  particle  with  given  x,  z,  and  v*  is  given  by 

uz^  <  [PylM- {qlMc)Ay{x,z)Y  -  [Pyl^- {qlMc)Ay{x,z  =  O)]^.  (C5) 

We  note  that  since  we  have  written  C5  without  making  explicit  use  of  the  Harris  model  of 
the  magnetic  field,  we  may  use  other  field  reversal  models  by  simply  using  the  relevant  model 
for  Ay  provided  that  the  approximations  regarding  x  and  Vx  are  still  valid. 
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Figure  1:  Modified  Harris  magnetic  field  geometry 


Figure  2:  X-line  magnetic  field  geometry 
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Figure  3.  Self-consistent  density,  current,  magnetic  field  and  pressure  tensor  profiles  calculated 
using  a  kappa  distribution  (k  =  4.5 )  with  a  bulk  drift  in  the  tailward  direction  of  =  0.1  u,* 
The  units  on  the  pressure  tensor  elements  are  nanodyne/cm^. 


Figure  4.  Plot  of  the  left  hand  side  of  equations  19  a-c  using  the  results  depicted  in  Figure  3. 
Since  all  three  curves  are  essentially  flat,  we  see  that  pressure  balance  is  maintained  to  a  high 
degree  of  accuracy. 
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Figure  5;  a)  Trajectory  of  an  ion  with  energy  H  =  0.7Al}9,^,  initial  pitch  angle  /3  =  100  , 
initial  phase  angle  (f)  =  0°  and  initial  starting  position  zjL  =  6.  b)  Current  profile  of  the 
particle  in  Fig.  3.a  as  calculated  from  the  standard  ridding  algorithm  with  the  top  grid  cell 

located  at  z/L  =  5. 


Figure  6:  Current  profile  for  the  particle  in  Fig.  5a  at  the  three  times:  a)r 
h)T=  400  and  c)T  =  600  The  top  grid  cell  is  taken  to  htz/L  =  8. 


Figure  7:  a)  Vyj^oL  as  a  function  of  z/L  for  the  particle  in  Fig.  5a.  The  particle  is  non 
midplane  crossing  on  the  segments  a-b  and  c-d,  and  we  use  guiding  center  coordinates.  The 
particle  is  midplane  crossing  on  the  segment  b-c  and  we  use  true  particle  coordinates.  Notes 
that  if  actual  particle  coordinates  are  used  on  the  segments  a-b  and  c-d,  there  is  an  excess  of 
positive  y-velocity  in  the  regions  5.3  <  z/L  <  6.7  and  -6.3  <  z/L  -  5  and  an  excess  of 
negative  y-current  in  the  regions  0.5  <  z/L  <  1.9  and  -2  <  z/|  <  0.3.  b)  y/L  as  a  function 
of  z/L  for  the  same  particle.  Note  the  VB  drift  in  the  regions  around  |z/L|  =  2  and  the  large 
y-velocity  at  |z/L|  =  1. 
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Figure  8;  Schematic  diagram  of  ion  and  electron  gyro  orbits  in  a  constant  magnetic  field 
B  =  Bqk.  Note  that  above  the  guiding  center  the  particles  have  only  positive  y-current, 
whereas  below  the  guiding  center  the  particles  have  only  negative  y-current 
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Figure  9:  Magnitude  of  the  positive  y- velocity  deposited  near  the  launching  point  of  the  particle 
as  a  function  of  the  pitch-angle  and  for  (a)  =  0.1,  (b)  =  0.2  and  (c)  b^  =  0.4.  The  solid 

line  is  the  approximate  expression  given  by  (22a)  and  (22b).  The  dots  are  found  by  pushing 
the  particles  through  the  magnetic  field.  The  dashed  vertical  lines  indicate  the  pitch-angle  that  * 
below  which  there  is  no  net  positive  velocity  left  at  the  launching  point. 
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Figure  10:  Plot  of  the  critical  pitch-angle  Pc  for  which  the  error  current  at  the  midplane  equals 
the  meandering  current  at  the  midplane  as  a  function  of  as  calculated  from  Eq.  24  (•).  The 
other  points  (■)  are  calculated  using  the  numerically  determined  values  of  {vy)p. 
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Figure  11.  Self-consistent  density,  current  and  magnetic  field  profiles  calculated  using  a 
Maxwellian  with  a  bulk  drift  in  the  tailward  direction  of  vd  =  0.1  Vth  and  a  high  energy 
isotropic  tail  which  comprises  10  percent  of  the  total  density. 
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Figure  12.  Self-consistent  density,  current  and  magnetic  field  profiles  calculated  using  a 
Maxwellian  source  distribution  of  particles  with  a  bulk  drift  in  the  tailward  direction  of  ^ 

vd  =  2  Vth- 
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Figure  13:  Approximate  particle  orbits  used  in  deriving  Eq.  (A.9-A.1 1) 


Figure  14:  Plot  showing  graphical  solution  to  equation  (A.5).  TVo  new  roots  to  the  equation 
occur  when  the  amplitude  of  the  sinusoidal  oscillation  fe”^tan(^fc)  sin(\/l  +  t)  increases 
sufficiently  to  cross  the  line  r.  The  curve  depicted  are  for  (a)  Pk  =  tan“^  {b^},  (b) 

Pk  =  tan~^i5irbJ2),  and  (c)  Pk  =  tan"^(97rfez/2). 
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Figure  15:  Effective  potential  for  the  z-motion  for  (Py/M—  ilnx)/{^oL)  =  (a)  1,  (b)  0,  and 
(c)  -1.  The  dashed  lines  correspond  to  orbits  of  different  energies  as  compared  to  effective 
potential  (c).  Particle  A  violates  the  inequality  in  (C4)  and  is  midplane  crossing.  Particle  C 
represents  a  solution  to  the  inequality  (C4)  and  is  non-midplane.  Particle  B  represents  the 
equality  condition  in  (C4)  and  is  at  the  transition  point  between  crossing  and  non-crossing 
motion. 
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